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ABSTRACT 

We analyze two multi-chord stellar occultations by Pluto observed on July 18th, 2012 and 
May 4th, 2013, and monitored respectively from five and six sites. They provide a total of fifteen 
light-curves, twelve of them being used for a simultaneous fit that uses a unique temperature 
profile, assuming a clear (no-haze) and pure atmosphere, but allowing for a possible pressure 
variation between the two dates. We find a solution that fits satisfactorily (i.e. within the noise 
level) all the twelve light-curves, providing atmospheric constraints between ~1,190 km (pressure 
~11 /ibar) and ~ 1,450 km (pressure ~0.1 fi bar) from Pluto’s center. Our main results are: (1) 
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the best-fitting temperature profile shows a stratosphere with strong positive gradient between 
1,190 km (at 36 K, 11 /ibar) and r = 1,215 km (6.0 /rbar), where a temperature maximum of 
110 K is reached; above it is a mesosphere with negative thermal gradient of -0.2 K km -1 up to 
~1,390 km (0.25 /ibar), where, the mesosphere connects itself to a more isothermal upper branch 
around 81 K; (2) the pressure shows a small (6%) but significant increase (6-tr level) between the 
two dates; (3) without a troposphere, Pluto’s radius is found to be Rp = 1,190±5 km. Allowing 
for a troposphere, Rp is constrained to lie between 1,168 and 1,195 km; (4) the currently measured 
CO abundance is too small to explain the mesospheric negative thermal gradient. Cooling by 
HCN is possible, but only if this species is largely saturated; Alternative explanations like zonal 
winds or vertical compositional variations of the atmosphere are unable to explain the observed 
mesospheric trend. 

Subject headings: planets and satellites: atmospheres, planets and satellites: physical evolution, 
methods: data analysis, methods: observational, techniques: photometric 


1. Introduction 


Stellar occultations are a very powerful tool to discover and study, among others, tenuous atmospheres 


around remote bodies. Pluto’s atmosphere was discovered using this technique (Brosch 1995 Elliot et al. 


1989 

Hubbard et al. 

1988), and its spectacular two-fold expansion between 1988 and 2003 was also revealed 

using stellar occultations (lElliot et al. 2003 

Sicardy et al. 2003 

. Other trans-neptunian objects were explored 


with this technique, and so far, none of them exhibited atmospheres at the 10 nbar pressure level (three 
orders of magnitude smaller than for Pluto). This includes Charon (Sicardy et al. 2006 ), Eris (Sicardy et al. 
2011), Makemake ( Ortiz et al.||2012 ) and Quaoar ( Braga-Ribas et ah||2013 ). 


All those bodies have sizes and surface gravities that are comparable to those of Pluto, within a factor 
of two. As such, the derived upper limits constrain the physical conditions necessary for the appearance and 
maintenance of atmospheres around a body with a given ice composition and heliocentric distance. 


Here we analyze results derived from two Pluto stellar occultations (18 July 2012 and 04 May 2013) 
that provide signal-to-noise ratios (SNR’s) that are among the best ever obtained during such events. They 
are furthermore combined with well-sampled multi-chord coverages, providing a good absolute radial scale 
for the atmosphere extension. 


We use the simplest possible model, assuming a spherically symmetric, clear (no-haze), pure N% atmo¬ 
sphere with constant temperature both horizontally and with time. Our model satisfactorily fits twelve of 
the selected light-curves, and provides accurate density, pressure and temperature profiles for radii between 
1,190 km (11 /ibar pressure level) and 1,450 km (~ 0.1 /ibar) from Pluto’s center, while also providing 
constraints on Pluto’s radius. 


As Pluto’s atmospheric pressure is dominated by the vapor equilibrium pressure at its surface, it is 
very sensitive to tiny changes of temperature and the available amount of exposed ice. This induces strong 


1 Partly based on observations made with the ESO camera NACO at the Very Large Telescope (Paranal), under programme 
ID’s 089.C-0314(C) and 291.C-5016. The prediction use observations made with the WFI camera at the 2.2 m Telescope, under 
programme ID’s 079.A-9202(A). 
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seasonal effects over the plutonian year (Hansen and Paige 1996) that can be monitored and analyzed through 
stellar occultations ( Young|2 013). In that context, our data reveal a small, but significant increase of pressure 
between 2012 and 2013, which can be used for constraining current Pluto seasonal models, see |01kin et al.| 
(2015) for a detailed analysis. 


Our results are obtained in the context of the forthcoming flyby of the dwarf planet by the NASA New 
Horizons spacecraft in July 2015. Consequently, they can be used as a basis of comparison with the New 
Horizons findings. 


2. The 2012 and 2013 Pluto stellar occultations 
2.1. Predictions 


From astrometric observations along Pluto’s path onto the sky plane between 2008-2015, performed at 
the ESO’s 2.2 m telescope, Assafin et al. (2010) made accurate predictions for stellar occultations involving 
the dwarf planet and its satellites. 


In this context, the two occultations analyzed here, one on 18 July 2012 and the other on 04 May 2013, 
stood out as promising events, owing to the magnitudes of the candidate stars and to the presence of several 
potential observing sites along the shadow’s path. 


Follow-up astrometric observations of the stars were carried out in order to improve the predictions. 
These observations were made with the 1.6 m (Perkin-Elmer) and 0.6 nr (Boiler & Chivens) telescopes, at 
Pico dos dias Observatory (OPD, IAU code 874), and they are done wherever possible within our access 
time. 


Moreover, 16 positive detections of other occultations by Pluto, that occur between 2005 and 2013, were 


used to improve Pluto’s ephemeris offset (see Benedetti-Rossi et al. (20141 for details). 


Days before the event, we carried out observations with Pluto and the occulted star present in the same 
field of view of our CCD’s in order to minimize systematic biases like catalog errors. 


2.2. Observations 

The 18 July 2012 Pluto occultation was observed near zenith from five sites in South America (Fig. [l]). 
The circumstances and technical details of the observations are provided in Table |T| The 04 May 2013 event 
was recorded from six sites, under similar conditions (Fig. [I]), providing ten light-curves (Table [2]). Various 
astrometric, photometric and physical parameters associated with each event are summarized in Table [3] 

Fig. m displays the reconstructed geometries of each event, showing the Plutocentric latitudes and 
altitudes probed by the primary stellar image at each site, see appendix for details. For Paranal (18 July 
2012), we plot for sake of illustration the trajectories of both primary and secondary stellar images. As 
commented later, the contribution of the secondary image is small but not negligible compared to that of 
the primary image near the shadow center. Note that the primary image probes the summer (resp. winter) 
Pluto’s hemisphere at ingress (resp. egress). 
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2.3. Photometry and calibration 


Classical bias, dark, fiat-field and sky subtraction provide the occultation light-curves displayed in 
Figs. [3j [3] and [5j In all cases, a reference star brighter than the target was used to correct for low-frequency 
transparency variations. 


As expected, the best SNR light-curve was obtained at Paranal on 18 July 2012, using the NAOS- 
CONICA0 (NACO) ( |Lenzen et al. 2003 Rousset et al. 2003) camera attached to the 8.2-m “Yepun” Very 


Large Telescope (VLT) of the European Southern Observatory (ESO), at a rate of 5 frames per second in 
H band. Moreover, this is the only data set for which we have an accurate photometric calibration, which 
allows us to subtract the contribution of Pluto and Charon from the occultation light-curve, see below. As 
such, the 18 July 2012 data provide the best constraints on Pluto’s atmospheric structure. However, the 04 
May 2013 light-curves have on average better SNR’s than those of 18 July 2012, as well as a better spatial 
sampling, thus providing better constraints on the absolute vertical scale of the atmosphere. 


Calibration images were taken with NACO some twenty minutes before the 2012 event. They show 
resolved images of Pluto, Charon and the star under excellent seeing conditions (Fig. [6]) . Digital coronagraphy 
(Assafin et al. 2008 2009) was used to remove the star contamination from Pluto and Charon images. 
Classical aperture photometry finally provided the Pluto + Charon flux relative the occulted star. This 
allows us to estimate the residual stellar fux in the deepest part of the 18 July 2012 occultation at Paranal, 
with a value that varied from 2.3 ± 0.8% to 1.8 ± 0.8% of its unocculted value in the central part of the 
occultation (Fig. [7]) . 


3. Modeling of Pluto’s atmosphere 


The general idea for modeling Pluto’s atmosphere is to use an iterative procedure, combining both direct 
ray-tracing and inversion approaches. We first invert our best signal-to-noise ratio light-curve to retrieve 
Pluto’s atmospheric density, pressure and temperature profiles (see Appendix and Vapillon et al. (1973)). 
The retrieved temperature profile is then used as a guide to generate, through direct ray-tracing, synthetic 
occultation light-curves that are simultaneously fitted to all the observed light-curves obtained at a given 
date. This pins down the location of Pluto’s shadow center relative to the occultations chords for both the 
2012 and 2013 events (Fig. [ 2 ]). Finally, the inversion of the best light-curve is performed again and the 
procedure is resumed. This iterative process eventually provides the accurate geometry of each event, as 
well as consistent density, pressure and temperature profiles that best fit all the occultation light-curves. 


Simplifying assumptions are made in our procedure (possible caveats are discussed later): (?) Pluto and 
its atmosphere are spherically symmetric, all quantities depend only on the radius r (defined as the distance 
to Pluto’s body center), ( ii ) the atmosphere is transparent (no haze present), and (ra) it is an ideal gas 
in hydrostatic equilibrium, in our case, a pure molecular nitrogen N 2 atmosphere, neglecting other minor 
species like methane. Moreover, ( iv ) we assume that T(r) is time-independent, i.e. the temperature profiles 
are the same in 2012 and 2013. Once T(r) is derived as detailed later, the density and pressure profiles n(r) 
and p(r) are derived from the hydrostatic and ideal gas equations (Eq. |A2| ), once a boundary condition is 
provided, i.e. the pressure at a given radius. ( v ) Although T(r ) is taken as time independent, the pressure 
is not. This is justified by the fact that the pressure is very sensitive to Pluto’s surface temperature through 


2 NAOS-CONICA is Nasmyth Adaptive Optics System (NAOS) and Near-Infrared Imager and Spectrograph (CONICA) 
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the vapor pressure equilibrium equation. For instance a 1-K temperature increase at the surface results in a 
two-fold increase of pressure or so (Fig. [8]). Thus, the pressure is a free parameter in our fits. More precisely, 
Eq. A2 requires a boundary condition, once T(r ) is fixed. So, we use as a free parameter the pressure p r at 
an arbitrary radius r. We choose r = 1, 275 km for an easier comparison with other works that provide the 
pressure at that level (see e.g. Olkin et al. (2015)). This level corresponds to a normalized stellar flux of 
« 0.45 in the shadow plane. Once £> 1,275 is given, the density and pressure profiles n(r) and p(r) are uniquely 
defined. 


We choose the 18 July 2012 occultation light-curve obtained at VLT/NACO to perform the first in¬ 
version. We use this particular light-curve because it has the highest SNR of all (Fig. [3]), and also because 
this is the only one for which we have a reliable measurement of the background contribution from Pluto 
and Charon (Fig. [7]), necessary to correctly invert any occultation light-curve. The successive steps of our 
procedure are as follows: 


(1) The inversion reveals a strong increase of temperature just above the surface (stratosphere), followed 
by a turning point where the temperature reaches a maximum (stratopause), then a region with a mild 
negative gradient (mesosphere), and finally an isothermal upper branch, see Appendix and Fig. |13| Using 
the prescriptions described by Eqs. Al we adjust the coefficients ci, ...c g controling the profile T(r) in order 
to best fit the inverted temperature profiles (see Table [4]) . 

(2) Keeping the profile T(r) fixed in shape, we simultaneously fit seven of the light-curves obtained 
on 04 May 2013. The free parameters of that fit are the two coordinates defining the shadow center, the 
pressure Pi ,275 at radius 1,275 km, and the value of rq, the deepest point that we consider in our profile. 
At this stage, when rq is varied, all the other radii rq, r 3 and rq defining T(r), see Eq. Al are changed by 
the same amount. In other words, the entire profile T(r) is vertically displaced by this amount. Thus, rq 
eventually fixes the absolute vertical scale of the atmospheric profile. Note that rq is not , a priori, the radius 
of the stratobase, nor Pluto’s surface radius. In practice, the choice of rq is made so that the stellar rays 
from the faint secondary image passing at rq have a contribution to the total flux that is negligible compared 
to the light-curve noise level. Thus, taking larger values of rq would create artificial discontinuities in the 
synthetic light-curve, while smaller values would require useless computation time. To find Pluto’s shadow 
center, we separate the fit along the direction of the star motion relative to Pluto from the fit perpendicular 
to that direction. This is because the fit along the star motion is essentially independent of the atmospheric 
model, and is generally more accurate than the fit perpendicular to that direction. 


Note that the 2013 light-curves have generally a better SNR than those of 2012 (excluding the VLT data 
set), because of a better distribution of the chords (Fig. [ 2 ]). Consequently, the 2013 occultation light-curves 
provide a better constraint for rq, or equivalently, for the absolute vertical scale of the atmospheric model, 
than those of 2012. 


(3) Fixing ?q to its value found in step (2), we turn back to the 18 July 2012 data set and simultaneously 
fit the five corresponding light-curves, varying Pluto’s shadow center and ^ 275 . 

The procedure is then resumed at point (1). It is a converging process that provides consistent solutions 
for the shape of the profile T(r), the absolute vertical scale for T(r), the centers of Pluto’s shadow for both 
events, and the two boundary conditions pi ,275 for the 18 July 2012 and 04 May 2013 events. This fitting 
procedure has a total of twelve free parameters: the nine coefficients C\, ...Cg, the two coordinates that define 
Pluto’s shadow center, and the pressure £> 1 , 275 - 

As commented before, the 18 July 2012 NACO light-curve is the only one for which the Pluto + Charon 








- 7 - 


contribution is measured (Fig. [7|. So, the stellar flux was normalized between that value and the full 
unocculted flux before starting the fit procedure. 


For the other light-curves, the inverse approach was used: the background Pluto + Charon flux was 
imposed by linearly adjusting the normalized, synthetic stellar flux to the actual occultation light-curve, 
through a least-square fit. As the residual stellar flux well inside the shadow is mainly controlled by the 
density scale-height of the deep stratosphere (Eq. B4), this means that the structure of that region is in 
fact dominated by the NACO, 18 July 2012 data. The other light-curves thus mainly serve to constrain the 
atmospheric structure above that level (mesosphere). 


4. General atmospheric structure 


The best fits of our synthetic light-curves to the data are shown in Figs. [3] and [4] For each light-curve, 
the residuals are displayed at the bottom of the corresponding panel. They show that a unique global model 
satisfactorily explains all the observations, with y 2 values per degree of freedom (y 2 iof , see Eq. A8) close to 
unity, except for the 18 July 2012 NACO data (Figs. [3] and [4]). In fact, due to the quality of this particular 
data set, the residuals are dominated by spikes associated with wave activity, as illustrated in Fig. [7] The 
wave activity, including the one observed in the NACO data, is discussed in details elsewhere, see |French| 


et al. (2015). 


The parameters of the best atmospheric model are listed in Table [4] Note that the only parameter 
that differs between 18 July 2012 and 04 May 2013 is the boundary condition, i.e. the pressure £> 1,275 at 
r = 1,275 km. Table [4] reveals a small ( 6 %) but significant ( 6 -a level) increase of pressure, from £> 1,275 = 
2.16 ± 0.02 /rbar in July 2012 to £> 1.275 = 2.30 ± 0.01 £tbar in May 2013, corresponding to a pressure increase 
rate of 7.5% per year. 


Based on various occultation data collected, Young (2013) and Olkin et al. (2015) report a general 
pressure increase of some 3.5-7.5% per year between 2006 and 2013, consistent with our result above. Note 
that our value of £> 1,275 = 2.30 ± 0.01 bar for May 2013 differs from Olkin et al. (2015)’s result (2.70 ± 0.2 
£tbar) by a barely significant 0.4±0.2 /rbar. Part of this difference could be due to the different methods used 
to derive those numbers, as Olkin et al. (2015) use an isothermal fit to the upper part of the light-curves, 


while we use a combination of mesosphere with negative thermal gradient and an upper isothermal branch 

(Fig.pl. 


Fig.[8]displays the density vs. radius, the temperature vs. pressure, the temperature vs. radius, and the 
temperature gradient of our best model. Also shown superimposed in that figure are the ingress and egress 
profiles retrieved from the inversions of the 18 July 2012 NACO light-curve, corresponding respectively to 
the summer and winter hemispheres, as far as the primary stellar image is concerned. Fig. [9] is a more 
detailed view of the bottom of the temperature and temperature gradient profiles, close to Pluto’s surface. 


The shaded areas in Fig. [8] and [ 9 ] indicate the l-cr error envelopes caused by (i) the photometric noise 
in the NACO light-curve, that mainly affects the upper parts of the profiles, and (ii) the uncertainty on the 
Pluto + Charon contribution to the total observed flux, that mainly affects the lower parts of the profiles. 
The methods to calculate these uncertainty domains are described in the Appendix. The temperature profiles 
are furthermore affected by another source of uncertainty, namely (in) the a priori unknown temperature 
boundary condition, inherent to the nature of Eq. |A2| (a first order differential equation). As examples, we 
show in Fig. [8] (gray lines in panels b, c and d) the profiles obtained by changing by ±5 K the nominal 





















boundary condition (T = 80.5 K at r = 1, 390 km) of the egress, inverted NACO temperature profile. 


Both the photometric noise and the unknowledge of the temperature boundary condition cause an 
exponential divergence of the uncertainty domain for T and dT/dr as r increases, with an e-folding distance 
equal to the density scale height H (Eq. B3 and Fig. [8]). Nevertheless, we note that if we have independent 
information on Pluto’s atmosphere, e.g. from theoretical models or forthcoming observations from the New 
Horizons mission, then we can constrain our temperature at rather high altitudes. For instance, at radius 
r = 1,450 km (pressure ~0.1 /ibar), the l-cr uncertainty on T caused by photometric noise is about ±2.5 K. 
Conversely, the two alternative solutions T{r) given as examples in Fig. [8] (the gray lines in panel (c)), using 
different boundary conditions, differ from each other by 30 K at that same radius. Consequently, they can 
be distinguished well above the noise level if we dispose of independent constraints on the thermal properties 
of the atmosphere at that radius. 


For instance, the warmer gray profile with strong positive temperature gradient in Fig.[8]can be discarded 


if we adopt current models which predict that UV heating is efficient only at much higher levels (Zhu et al. 


2014). The same is true for the cooler gray temperature profile in panel (c) of Fig. [8j it shows too strong a 


negative gradient in the 1,400-1,450 km range, considering that atmospheric escape may cause temperatures 
as low as ~ 60 K, but only much higher in the atmosphere (Ibid.). 


Figs. [H] reveals three regions in our thermal profile, from bottom to top: a stratosphere with strong 
positive gradient that starts around 1,190 km with temperature near 36 K and pressure 11 /ibar, and 
reaches a maximum temperature of 110 K at the stratopause (near r = 1,215 km, 6.0 /ibar). Then follows 
a mesosphere with mild negative thermal gradient of -0.2 K km -1 up to the mesopause (r ~ 1,390 km, 
0.25 /ibar), where, it connects itself to a more isothermal upper part around 81 K. These regions are now 
described in detail. 


5. Stratosphere 

As explained in the Appendix, the residual stellar flux in the mid-part of the occultation is proportional 
to the local density scale-height H, which is itself related to the strong stratospheric temperature gradient 
(Eq.|B4j. 

It is important to note that at closest approach to Pluto’s shadow center on 18 July 2012, our model 
predicts that the secondary image observed at Paranal contributes by 20% to the total, primary ± secondary 
stellar residual flux (Fig. [7]). This is not negligible and explains why we have to extend our ray tracing model 
below the deepest radius obtained for the inverted temperature profiles (red and blue lines in Fig.[9|. In fact, 
the inversion procedure assumes that there is only one (primary) stellar image contributing to the flux at 
any moment, while the direct ray tracing procedure does account for the presence of the two images. When 
the secondary image appears and disappears (at the extremities of the orange trajectory shown in Fig. 
reaches the radius rq = 1,190.4 km (Table [4]). Its appearance and disappearance cause small discontinuities 
in the synthetic flux, but they are too small to be distinguished from the noise (Fig. [7|. 

Due to the uncertainty on the Pluto ± Charon flux contribution, the deepest point of our model is 
determined to be at 1,190 ± 5 km (Fig. [9]). At that point, nitrogen reaches its saturation vapor pressure 
(Fig-i. and thus condenses in principle into ice, i.e. reaches Pluto’s surface. In that context, we obtain a 
solution with a clear nitrogen atmosphere and a Pluto radius of 1,190 ± 5 km which consistently explains 
all our observations, accounting for the presence of both the primary and secondary images. 
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Other models are possible, though. Based a more incomplete and lower quality data set than used 


here, Lellouch et al. (2009) conclude that the nitrogen condensation level occurs somewhere in the range 


1,187-1197 km, consistent with the present work. However, a shallow adiabatic troposphere with dry or wet 
nitrogen (or methane) may exist below 1,190 km. Nevertheless, there is little freedom for such tropospheric 
models because (i) they tend to create caustics in the light-curves that are not observed and (ii) they provide 
a cold methane column density that would be detected by other means. More precisely, using spectral data, 


Lellouch et al. (2009) find possible tropospheric solutions in a narrow region of the parameter space, with 


depths that cannot exceed 17 km. Similarly, combining again constraints from spectra with a preliminary 


analysis of the occultation data presented in this work, Lellouch et al. (2015a) concluded that Pluto’s radius 


should be between 1,180-1,188 km, some 2-8 km below the condensation radius 1,190 km derived above. 
This said, we assume here that the atmosphere is haze-free, a subject of debate since the discovery of 


Pluto’s atmosphere. Analyzing a high SNR occultation observed in 2006, Young et al. (2008) conclude that 


a haze-only explanation for the light-curve is extremely unlikely. In fact, the clear atmosphere model implies 
a temperature profile that naturally connects the maximum temperature of ~110 K near 1,215 km to the 
surface at average temperature of ~50 K (Lellouch et al. 2000 2013), see Fig. [8]. 

Other constraints come from a central flash observation during a stellar occultation in July 2007. From 


that event, Olkin et al. (20141 conclude that the flash is consistent with a transparent atmosphere with 
temperature gradient of 5 K km” 1 at 1,196 km, fully consistent with our results (Fig. Oj). Olkin et al. (2014) 


exclude in particular a haze-only model to explain the central flash, although combinations of thermal 


gradient and haze mechanism are possible. In the same vein, Gulbis et al. (2015) use a wavelength-resolved 


occultation on 2011 to constrain the presence of hazes in Pluto’s atmosphere. Although haze models do 
improve the fit residuals, a clear atmosphere with a steep thermal gradient at the bottom is also consistent 
with the observations. 

Finally, we note that the residual stellar flux exhibits a significant decrease in the bottom of the light- 
curve, from 2.3% to 1.8% of its unocculted value, in the central part of the occultation as observed from 
Paranal on 18 July 2012 (Fig. [7]). This behavior was already pointed out by Sicardy et al. (2003), based on 
another hight SNR occultation observed in August 2002. In both cases, the residual stellar flux decreased as 
the primary stellar image scanned first the summer, permanently lit northern lower atmosphere, and then 
the winter, low insolation region (Fig. [2]). This point is discussed in the last section. 


6. Mesospheric negative temperature gradient 


Above the stratopause (r ~ 1, 215 km), the temperature profile exhibits a negative temperature gradient 
up to r ~ 1,390 km, with an average value of dT/dr ~ —0.2 K km” 1 . In this ~ 170 km radius interval, 
the temperature decreases by some 30 K. This mesospheric gradient is little affected by the choice of the 
temperature boundary condition, see Fig. |8j panel (d). While the photometric noise and the boundary 
condition problem induce rapidly diverging solutions for dT/dr above ~ 1,400 km, the thermal gradient 
between ~ 1,250 km and ~ 1,360 km is robustly constrained around -0.2 K km” 1 , with a typical fine-scale 
scatter of ±0.05 K km” 1 that is dominated by Pluto’s wave activity, and not by the photometric noise. 
In this interval, the thermal gradient remains smaller (in absolute value) than the dry adiabatic lapse rate 


- 9/cp (Fig. §> where g is the acceleration of gravity (Eq. |A3|) and c p = 1.04 x 10 3 J K 1 kg 1 is the specific 


heat at constant pressure for N^. Thus, the mesosphere remains convectively stable. 


Note that in principle we may choose an extreme temperature boundary condition that provides an 
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isothermal mesosphere, i.e. a thermal profile that is much warmer than the warmer gray profile shown in 
panel (c) of Fig. [ 8 ] As commented earlier, however, this replaces one problem by another one, namely that 
the upper part of our profile is too warm, with seemingly no plausible physical explanation. 


The negative mesospheric thermal gradient is further confirmed by the inversions of our best SNR 
light-curves obtained in July 2012 and May 2013, see Fig. [To] This eliminates random, low frequency sky- 
transparency variations that may have corrupted the light-curves. Moreover, such gradients have also been 
reported in previous, independent works. For instance Young et al. |2008) derive and discuss a dT/dr = 
—0.086 ± 0.033 K km -1 gradient at r — 1,275 km from the 12 June 2006 occultation, while Elliot et al. 
(2007) give —0.17 ± 0.05 K km -1 for the same occultation and at the same radius. Gulbis et al. (2015) 
report a gradient of —0.23 ± 0.05 K km -1 in the 1,310-1,450 km region from the 23 June 2011 occultation, 
consistent with Person et al. (2013) for that event. Finally, Bosh et al. (2015) derive values of —0.17±0.03 K 
km ^ 1 and —0.24 ± 0.01 K km -1 around 1,280-1,300 km, for occultations observed on 09 September 2012 
and 04 May 2013, respectively. 


The origin of this thermal gradient is still debated. Two classes of possible explanations can be proposed: 
(1) the presence of cooling minor species and (2) yet unmodeled physical mechanisms. They are now examined 
in detail. 


6.1. Possible cooling by CO or HCN 

Radiative-conductive models of Pluto’s atmosphere have been developed initially by [Yelle and Lunine 


(1989); Hubbard et al. (1988); Lellouch (19941; Lellouch et al. (2015a) to explain the recently-discovered 


gross characteristics of Pluto’s atmosphere: a large lower atmosphere temperature gradient, and a warmer 
(~100 K) mesosphere. These studies used a simplified description of the heating/cooling properties of Pluto’s 


atmosphere proposed by Yelle and Lunine (19891, with heating in the methane 3.3 /im band and cooling 


in its 7.6 /im band, both occurring in non-LTE conditions. Lellouch (19941 first suggested that additional 


cooling due to LTE CO emission rotational lines was important, based on an estimated abundance of CO in 
Pluto’s atmosphere (10 - 4 -10 -3 ). 


These studies were updated with the much more extensive model of Strobel et al. (1996). Notably these 


authors improved the treatment of solar heating in the CH 4 near-infrared bands by considering the effects 
of opacity and vibrational (V-V and V-T) energy transfer, and showed the need to include heating from the 
2.3 /xm band system in addition to the 3.3 /im bands. 

As the composition of Pluto’s atmosphere, as well as surface (pressure, radius) conditions, were largely 
unconstrained at that time, Strobel et al. (1996J explored diverse combinations of surface pressure and 
methane mixing ratios (including non-uniform ones), including also the effect of CO cooling. In general 
these models were reasonably successful at explaining large near-surface temperature gradients, though ( i ) 
fitting 10-20 K km -1 gradients near the surface required pushing the models to their limits, e.g. a 3.6% 
CH 4 mixing ratio confined to the first scale height near the surface and a 3 /ibar surface (or tropopause) 
pressure; (ii) models tended to overestimate the upper atmosphere temperature (~ 130 K instead of 100 K). 


A general feature of the Strobel et al. (1996) models was their prediction of a mostly isothermal atmosphere at 
pressures less than ~ 2 /ibar, though some models exhibited a moderate (0 to -0.1 K km -1 negative gradient 
at 1-2 /ibar. As the direct detection of N 2 in Pluto’s atmosphere is still missing, they also considered a CO- 
dominated atmosphere case (e.g. 97% CO + 3% CH 4 ). This case led, through enhanced CO cooling, to much 
larger negative temperature gradients in the sub-microbar region and an upper atmosphere temperature of 
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about 55 K. 


The avaibility of new, quantitative, observational constraints on the composition (CH 4 ~ 0.5%, CO ~ 
0.05%) and near-surface structure (surface radius and pressure, tropospheric depth) of Pluto’s atmosphere 
from near-IR observations ( Lellouch et al.|2009 2011 1 prompted a revival of the Strobel et al. (1996) models 
( Zalucha et al.|[2011a|b Zhu et ah|2014 ). 


Model updates included new estimates of the vibrational energy transfer based on recent laboratory 
measurements of collisional relaxation rates (Siddles et al. 1994 Boursier et al.||2003 ), as well as the intro¬ 
duction of a scheme parameterizing the processes of eddy mixing and convection. With the updated model, 


Zalucha et al. (2011a I explored the effect of parameter space (CH 4 and CO mixing ratios, surface pressure 


and radius) allowed by the recent observations, assuming uniform vertical mixing of CH4 and CO (which 
was recently demonstrated to be the case for CH 4 in the first 25 km of the atmosphere, Lellouch et al.| 


(2015a)). Radiative-convective calculations were then coupled to a model generating synthetic occultation 


lightcurves for direct comparison to observations. The study was extended by Zalucha et al. (2011b) to 
include a putative troposphere. 


earlier 

Strobel et al. 

(1996) models. The stratopause temperature is still somewhat too high (120-125 K) 

near 1,215 km radius 

m Zalucha et al. 

( 2011 a 

. These models generally show only weak negative temperature 

gradients above this level, typically a ~5 K decrease over a 300 km range for a CO mixing ratio of 5 x 10 4 , 


or mild ~10 K decrease due to atmospheric escape (Zhu et al. J2014J. This is in disagreement with the profile 
reported in the present study, which exhibits a typical 30 K decrease between 1,215-1,390 km and a gradient 
of -0.2 K km -1 , as discussed earlier. Rather, the profile we derive is remarkably similar to that calculated 
by Zalucha et al. (2011a) (their Fig. 8 ) for the case of 40-times enhanced CO mixing ratio (200 x 10 -4 ) This 


scenario, however, is at odds with the direct measurement of the CO abundance (Lellouch et al. 2011). This 
suggests that an additional cooling source is at work. 

Through radiation in its intense rotational lines, HCN is the major cooling agent in Titan’s upper 
atmosphere, where its mixing ratio is typically 2 x 10 -4 at 1,100 km (Vuitton et al. 2007), and where it 


equilibrates the solar UV heating rates (Yelle 19911. HCN has not been detected yet in Pluto’s atmosphere, 
but its presence is expected from coupled photochemistry in a N 2 -CH 4 atmosphere. A complete re-assessment 
of the Pluto models is however beyond the scope of the present study. Here we only re-calculate CO cooling 
rates, and also examine the case of HCN cooling. Photochemical models predictions lead to rather diverse 
mixing ratios of HCN (10 - 8 -10 -3 of N 2 (Summers et al. 1997 Lara et al. 1997 Krasnopolsky and Cruikshank 


1999), where the difference largely seems to come from the fact that the more “optimistic” models have not 


accounted for the fact that under cold (< 100 K) temperatures, atmospheric condensation of HCN should 
occur. Here, we nominally consider cases in which the HCN abundance is limited by the saturation law 


(Fray and Schmitt 2009), but we also run a case with uniformely mixed HCN, as supersaturation may be 


possible in a clear, tenuous atmosphere as Pluto’s. Cooling rates at radius r are calculated from the following 
equation (e.g. for CO): 

Rco{r) = 4:TrN C o(r) J B v {T(r))k v E 2 {T v )dv , ( 1 ) 

where Nqo is the local CO number density, T is temperature, k u and r„ are the absorption coefficients and 
zenithal opacity at frequency v, and E 2 (t) is the second exponential integral. The integral runs over the 
entire mm/submm range (0-200 cm -1 ), and unlike in Strobel et al. (1996), We include all isotopic variants 
of CO and HCN, i.e. lines of CO, 13CO, C180, HCN, H13CN and HC15N are taken into account when 
calculating the opacities. Moreover, absorption coefficients are calculated using a Voigt profile, instead of 
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the Doppler approximation. Both aspects lead to a minor but non-negligible increase in the cooling rate at 
low altitudes. Calculations of the cooling rates are performed for the thermal profile inferred in this work. 

Results are shown in Fig. |TT] for a series of assumed CO and HCN profiles. The CO mixing ratio 
QCO = 5 x 10 -4 curves show the “nominal” CO cooling. Although Zalucha et al. (2011a b) do not show their 
cooling rates, our calculation for CO can be compared to Figure 5a of Zhu et al. (2014), showing reasonable 


agreement. Increasing the qco by a factor of 40 leads to an increase of the cooling rate, although in much 
lower proportion due to opacity effects. 

Other curves in Fig.[lI]show the cooling due to HCN for different assumed HCN mixing ratios in the non- 
saturated part of the atmosphere. For the temperature profile derived in this work, the lower temperatures 
above 1270 km radius, severely restrict the amount of gaseous HCN if saturation of HCN is accounted for. 
In fact, HCN appears to be saturated everywhere in the atmosphere, except possibly over the 1210-1270 km 
range, where the condensation law allows 10 _7 -1CD 6 HCN mixing ratios. There, the HCN cooling rate may 


slightly exceed the nominal CO cooling rate (pink vs. red curves in the left panel of Fig. 11). However, 
for the HCN cooling rates to approach the “enhanced” CO cooling rates necessary to explain our negative 
mesospheric temperature gradient (i.e. those for qco = 200 x 10 -4 , as considered by Zalucha et al. (2011a)), 


one must assume that HCN is not limited by saturation. Specifically, the blue curve in Fig. El shows that a 
uniform HCN mixing ratio of ~ 5 x 10~ 5 is required. 

Although a full re-assessment of the radiative-models should be undertaken at this point, we conclude 
from this exercise that there is no obvious “culprit” for a -0.2 K km -1 temperature gradient above the 
radius ~ 1220 km. According to the calculations of Zalucha et al. (2011a), CO in amounts consistent with 


the direct observations of Lellouch et al. (2011) provide unsufficient cooling. We show here that HCN could 


be an alternative efficient cooling agent, but only if its mixing ration vastly exceeds expectations from the 
condensation law. Direct measurements/upper limits of HCN from ALMA and perhaps from New Horizons 
UV spectrometer (ALICE) will bring new light on this issuej^] 


6.2. Alternative explanations 


Coming back a step, the primary result derived from a stellar occultation light-curve is the refractivity 
profile z'(r), from which a density profile n(r) = v(r)/K is obtained, assuming a given gas composition 

(Eq.|A§. 

A first idea is to envisage that hazes are present in the mesosphere. Those hazes would absorb part of 
the stellar flux, in such a way that a basically isothermal mesosphere is thought to host a negative thermal 
gradient just because of the clear-atmosphere assumption. To test that hypothesis, we have generated 
synthetic light-curves, forcing the mesosphere to be isothermal at Xi so = 95.5 K above the stratopause (we 
have also tested other values of T iso between 85 and 110 K, with the same conclusions). Fig. [3] shows the 
resulting residuals for the NACO 18 July 2012 light-curve (labelled “iso.” in that figure). They depart from 
zero well above the noise level, and we can rule out photometric fluctuations caused by absorbing haze layers, 
since the residuals have both positive and negative values. 


This said, two assumptions may be wrong in Eq. A2 (1) the atmosphere may be not composed of pure 
nitrogen N 2 , so that the nitrogen molecular mass must be replaced by a new value //, and (2) hypothetical 


3 The detection of HCN in Pluto’s atmosphere, using ALMA, was announced by 


Lellouch et al. 


(2015b I on July 30, 2015. 


























- 13 - 


zonal winds may create a centrifugal acceleration, so that the acceleration of gravity g must be replaced by 
a term g' that includes supplementary terms. 


In fact, we can use Eq. A2 in a reversed way. More precisely, the refractivity profile v[r) is actually 


an imposed observable (since it is directly derived from the occultation light-curve), while we may use a 
prescribed temperature profile T'(r ), for instance taken from a theoretical model. With this approach in 
mind, Eq. | A2| can be re-written: 

pig' _ kT' d\og(vT') 

M3 M3 


dr 


( 2 ) 


To obtain this equation, we have used v[r ) = K ■ n(r), where K is the molecular refractivity (Eq. A4). We 
assume here that the atmospheric composition varies slowly with radius (i.e. K suffers small variations over 
one scale-height), so that we neglect dK/dr , and finally provide d{\og(v)\/dr ~ d[\og(n)} / dr. 

Thus, the ratio pig'/pg is the factor by which the molecular mass and/or the acceleration of gravity g 
must be multiplied in order to retrieve a prescribed temperature profile T'(r), given an observed occultation 
light-curve. 


In Fig. 12 we consider an example where the prescribed temperature profile T'(r ) exhibits a decrease 
of only 10 K between the stratopause and the mesopause. This is typical of what can be obtained by the 
combined effects of CO cooling ( Zalucha et al.pOlla ) and/or a atmospheric escape ( Zhu et aT7||2014 ). The 
right panel of Fig. 12 shows the resulting profile for pig'/pg, restricting ourselves to the mesospheric region. 


We first assume here that the atmosphere is composed of pure nitrogen, so that p' = p 1 and the ratio 
p'g' /pg = g'/g is only caused by variations of g' . In the presence of a zonal wind with velocity v, the 
centrifugal acceleration provides g' = g — v 2 /r, and from g = GM/r 2 , a zonal wind of: 


v = 


GM 


11 - 


9' 


840* 1 — — m s \ 


(3) 


where we have used the value of GM in Table [3] and r ~ 1,250 km. With the example above, the factor 
g'/g' reaches a minimum value of about 0.95, yielding v ~ 190 m s — . This is close to supersonic, as the 
speed of sound for nitrogen N 2 at 100 K is about 200 ms -1 . In fact, current General Circulation Models 
(GCM’s) for Pluto predict zonal winds of less than 10 m s _1 at the altitudes considered here (Vangvichith 
2013;|Zalucha and M ichaels| 2013). Moreover, we see that above r ~ 1, 300 km, the ratio g'/g becomes larger 
than unity with the example considered here, which is impossible from Eq. [3} Other prescribed profiles T'(r ) 
could be imagined to avoid this problem by displacing the p'g'/pg profile to the left in Fig. 12 (providing 
smaller values of g'/g), but this would imply even more unrealistic, high zonal winds. 

Considering that g'/g ~ 1 from the discussion above, the p'g'/pg profile would represent variations of 


molecular of the atmospheric molecular weight, pi/p. In the example of Fig. 12 the molecular weight of the 
atmosphere has to be inferior to that of molecular nitrogen, p , to mimic the effect of a negative temperature 
gradient. This could be caused by the presence of a lighter gas, for instance neon, which has a molecular 
weight ptie ~ 0.72 p. That species has a relatively large solar abundance ( N e /N ~ 1.5) and is not condensed 


at Pluto’s atmospheres temperatures. The minimum value p!/p = 0.95 near 1,230 km (Fig. 12) would then 
require a local neon abundance of about 82%. However, and as before, the ratio p'/p would be larger than 
unity above 1,300 km, requiring that another, heavier, gas (e.g. argon) takes over above 1,300 km and drives 
the molecular mass upwards. Such model is clearly unrealistic though, because mass separation would result 
in a depletion, not enrichment, of the heavier species in the upper atmosphere. 
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7. Upper atmosphere 

Above r ~ 1,400 km, the results of Fig. [8] shows that a change of thermal gradient may occur, with a 
more isothermal upper branch just above the mesosphere. However, the lack of independent constraints on 
the temperature at that level prevents an unambiguous choice for a particular solution for T(r). In addition, 
the rapidly increasing contribution of the noise makes impossible any estimation of the thermal gradient 
above 1,450 km (Fig. [8j panel (d)). This said, the lack of obvious mechanisms to drastically warm up or 
cool down the atmosphere just above 1,400 km suggests (but by no means proves) an isothermal branch 
between 1,400 and 1,450 km. Under this hypothesis, we estimate a 1-cr uncertainty domain of 81 ± 6 K for 
the temperature of this isothermal branch. (This interval corresponds to an increase of Ay 2 = +1 of the y 2 
function with respect to the best, minimum value y 2 lin ). 


8. Discussion and conclusions 


We have analyzed among the best light-curves ever obtained during stellar occultations by Pluto. Combi¬ 
nation of well-sampled occultation chords (Fig. i and high SNR data (Figs.|3p|) have allowed us to constrain 
the density, temperature and thermal gradient profiles of Pluto’s atmosphere between radii r ~ 1,190 km 
(pressure p ~ 11 fibax) and r ~ 1,450 km (pressure p ~ 0.1 pbar). Our main results are listed below. 

Global Pluto’s atmospheric model. We find that a unique thermal model can fit satisfactorily twelve 
light-curves observed in 2012 and 2013 (Figs. [3]and|4|, assuming a spherically symmetric and clear (no 
haze) atmoshere. The parameters defining our best model are listed in Table [4] (see also Fig. 13), and the 


various resulting profiles (density, temperature and thermal gradient) are displayed in Figs.pJO and 10 The 
absolute vertical scale of our global model has an internal accuracy of about ±1 km (Table HI). However, this 
error is amplified to ±5 km at the bottom of the profiles (Fig. [9]), because of the uncertainty on the residual 
stellar flux (Fig. [T]) in the central part of the occultation observed by NACO on 18 July 2012. 

We quantify in this work the propagation of the photometric noise into the density, temperature and 
thermal gradient profiles (Eqs. B3 and Fig. [8]). The key parameter that governs the noise propagation is 
the radius ro in the amosphere at which the stellar drop caused by differential refraction is equal to the flux 
standard deviation. The radius ro can be estimated from Eq. |B2[ which includes all the quantities at work 
in a stellar occultation: photometric noise, molecular refractivity, atmospheric scale-height and radius, and 
distance to the body. For the NACO light-curve, we find ro = 1,565 km, corresponding to a pressure level 
of about 14 nbar. 

Although a satisfactory fit to all the data used here is provided by a unique model, there are two slight, 
but significant departures from this global model. They are now discussed in turn. 

Pressure increase between 2012 and 2013. In the frame of our model (i.e. assuming a constant tempera¬ 
ture profile), we detect a significant 6% pressure increase (at the 6-cr level) during the ~9.5 months separating 
the two events under study. This means that Pluto’s atmosphere was still expanding at that time, confirming 


the work of Olkin et al. (2015), which compiles and analyzes pressure measurements between 1988 and 2013. 


Ingress/egress asymmetry of lower temperature profiles. Fig.[7]shows that the stellar flux decreased from 
2.3% to 1.8% of its unocculted value during the central part of the 18 July 2012 occultation, as observed by 
NACO from Paranal. This corresponds to the primary stellar image first scanning the summer, permanently 
lit Pluto northern hemisphere, and then the winter low insolation southern hemisphere (Fig. [2]). This 
confirms a similar trend pointed out by Sicardy et al. (2003) during another high SNR stellar occultation 
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recorded in August 2002. These authors interpreted this result as a surface boundary layer effect, where 
the lowermost scale-height adjusts itself to the surface temperature variegations, which might explain the 
behavior displayed in Fig. [9] 

Another interpretation of this trend is the gradual entrance of the primary stellar image into an absorbing 
haze layer near Pluto’s evening limb, a hypothesis that can be tested during the New Horizons flyby in July 
2015. 


Pluto’s radius and density. The extrapolation of our temperature profiles to the nitrogen saturation line 
implies that nitrogen may condense at a Pluto’s radius of Rp = 1,190 ± 5 km. However, the few kilometers 
above Pluto’s surface remain “terra incognita” as far as stellar occultations are concerned. In particular, 
the temperature gradients shown in Fig. [9] may deviate from the simple extrapolation used here, especially 
if haze layers affect the retrieved temperature profiles. Although difficult to envisage because of the strong 
caustics that they cause, a troposphere below 1,190 ± 5 km cannot be excluded. Combining high-resolution 
spectroscopic observations of gaseous methane, combined with constraints from an occultation observed in 


2002, Lcllouch et al. (2009|) conclude that the troposphere depth cannot exceed about 17 km. Consequently 


(and assuming that the temperature of the deep atmosphere did not change significantly since 2002), our 
observations constrain Pluto’s radius to lie in the range 1,168-1,195 km. More recently, combining constraints 


from spectra and a preliminary analysis of the occultation data presented in this work, Lellouch et al. (2015a) 


concluded that Pluto should have a radius between 1,180-1,188 km, some 2-8 km below the condensation 
radius of 1,190 km that we derive above. 


From a Pluto’s mass of Mp = 1.304 ± 0.006x 10 2 " kg (Tholen et al. 2008), we derive a density pp = 


(1.802 ± 0.007)(i?p/l, 200 km) -3 g cm -3 . Our estimation Rp = 1,190 ±5 km thus implies pp = 1.85 ±0.02 
g cm -3 in the absence of troposphere, and a range pp = 1.83— 1.95 g cm -3 if a troposphere is allowed. This 
is larger, but not by much, than Charon’s density, pc = 1.63 ± 0.05 g cm -3 (Ibid.). 

The mesospheric negative thermal gradient. Pluto’s stratopause occurs near 1,215 km (pressure p = 6.0 
/rbar), with a maximum temperature of 110 ± 1 K, where the error bar applies to the best inverted profile 
(NACO 18 July 2012), and stems from the uncertainty on the Pluto + Charon flux contribution (Fig. [8]). 

Above the stratopause, and up to about 1,390 km, our best 2012 and 2013 occultation light-curves yield 
inverted temperature profiles with a negative thermal gradient close to -0.2 K km -1 which amounts to a 
total decrease of 30 K for the temperature between 1,215 and 1,390 km (Figs. M 

Explaining this negative gradient by CO cooling requires a mixing ratio (200 x ICC 4 ) that is too high by 


a factor of 40 compared to current measurements (Lcllouch et al. 2011). Cooling by HCN is also discussed 


in this paper. It appears to be a possible alternative solution, but only if it remains largely supersaturated 
in the mesosphere. 


Changing the temperature boundary condition may suppress the negative gradient, but at the expense 
of creating a warm, unexplained thermal profile above 1,350 km. We have investigated more exotic solutions, 
like zonal winds or compositional variations that would “unbend” the retrieved temperature profiles, allowing 
a more isotliemal mesosphere. However, no realistic models could be built upon those alternative assump¬ 
tions. Again, the New Horizons flyby will provide constraints on the temperature boundary conditions and 
atmospheric composition that will be used to discriminate between the various solutions decribed here. 
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Fig. 1. The post-occultation, reconstructed paths of Pluto’s shadow for the two events studied here. The 
red dots indicate the shadow center every minute and the arrows show the direction of motion. The green 
dots mark the sites where data were obtained. The black solid lines correspond to the half-light stellar level, 
while the dotted lines correspond to the 1% stellar drop, thus marking the practical region of detectability 
of the occultations. Left - The 18 July 2012 event. The first red dot at right is at 04h 09m UT, the last one 
at left corresponds to 04h 18m UT. Right -The 4 May 2013 event. The first red dot at right is at 08h 12m 
UT, the last one at left corresponds to 08h 33m UT. 

























- 23 - 


18 July 2012 



Huancayo 


SP de Atacama 
Paranal 


Cerro Burek 
S Martina 


04 May 2013 



SP de Atacama 
Paranal 
Pico do Dias 


La Silla 
Cerro Tololo 
Cerro Burek 


Fig. 2.— Left - The trajectories of the primary stellar images relative to Pluto, as seen from the five stations 
used on 18 July 2012, see Table [l] The black arrow shows the general direction of stellar motion. Here, 
Pluto’s has an assumed radius Tip = 1,190 km (see text), and its center is indicated by the cross symbol. 
The gray arrow inside the disk indicates the direction of rotation. In the case of Paranal, we have plotted 
the path of the primary image in green, and the associated path of the secondary image in orange (see also 
Fig-0- The green and orange arrows show the corresponding local stellar motion along Pluto’s limb. Note 
that the two images move in opposite directions. The black star symbol shows the star position as seen 
from Paranal at a given, arbitrary moment, while the green and orange star symbols indicate the associated 
primary and secondary images at that time, respectively. Note that the three star symbols and the cross are 
aligned. Right - The same as left panel for the 04 May 2013 occultation, with only the paths of the primary 
stellar images plotted. In both panels, the summer, permanently lit Pluto’s hemisphere is at right, and the 
low insolation winter limb is at left. 
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Table 1: Circumstances of the 18 July 2012 Pluto occultation 


Site 

Lat. (dams) 
Lon. (dams) 
alt (m) 

Telescope 

Instrument/filter 

Exp. time/cycle (s) 1 

Observers 

Observatory UC 
(Santa Martina) 

33:16:09.0 S 

70:32:04.0 W 
1,450 

0.4 nr 

CCD/clear 

1 .0/1.0 

R. Leiva Espinoza 

Cerro Burek 

31:47:12.4 S 

69:18:24.5 W 
2,591 

ASH 2 0.45 m 
SBIG-STL11000/clear 

13.0/15.7 

N. Morales 

Paranal 

24:37:31.0 S 

70:24:08.0 W 
2,635 

VLT Yepun 8.2 m 
NACO/H 

0 .2/0.2 

J. Girard 

San Pedro 

de Atacama 

22:57:12.3 S 

68:10:47.6 W 
2,397 

ASH2 0.4 m 
SBIG-STL11000/clear 

13.0/15.44 

N. Morales 

Huancayo 

12:02:32.2 S 

75:19:14.7 W 
3,344 

0.2 nr 

CCD/clear 

10.24/10.24 3 
5.12/5.12 3 

E. Meza 


Notes. 1 Cycle is defined as the exposure time plus the readout time also known as dead time. Observations 
with the same exposure time and cycle have no dead time. 

2 ASH - Astrograph for the Southern Hemisphere. 

3 Exposure time was changed at 04:11:46 UT 
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Table 2: Circumstances of the 04 May 2013 Pluto occultation 


Site 

Lat. (dams) 

Lon. (dams) 
alt meters 

Telescope 

Instrument/filter 

Exp. time/cycle (s) 

Observers 

Cerro Burek 

31:47:12.4 S 

69:18:24.5 W 
2,591 m 

ASH 0.45m 

SBIG-STL11000/clear 

6/8 

J.L. Ortiz 

N. Morales 

CASLEO 

31:47:55.6 S 

Jorge Sahade 2.15m 

5/6.8 

R. Gil-Hutton 

(Leoncito) 

69:17:44.9 W 
2,492 m 

CCD/R 


C. Lopez-Sisterna 

Cerro Tololo 

30:10:03.4 S 

70:48:19.0 W 
2,207 

PROMPT 1 0.4m 

PI, P3, P4, P5 

CCD/clear 

5/8 

P3 offset 2 s 

P4 offset 4 s 

P5 offset 6 s 

J. Pollock 

La Silla 

29 15 21.276 S 

70 44 20.184 W 
2,336 

Danish 1.54m 

Lucky Imager/Z (A > 650nm 
CCD/iXon response) 

0 .1 /0.1 

several 
interruptions 
due to image 
cube writing 

L. Mancini 

La Silla 

29 15 16.59 S 

70 44 21.82 W 
2,315 

TRAPPIST 2 0.6m 

CCD/clear 

4.5/6 

E. Jehin, 

A. Decock, M. Gillon 
C. Opitom 

Pico dos Dias 

22 32 07.8 S 

45 34 57.7 W 
1,811 

B&C 3 0.6m 

CCD/I 

5/5.40 

M. Assafin, 

A. Ramos-Gomes Jr 

Ponta Grossa 

25 05 24.00 S 

50 09 36.00 W 

909 

Meade 16 0.4m 

CCD/clear 

5 

Tecnichal 

Problems 

M. Emilio 

Cerro Paranal 

24:37:31.0 S 

70:24:08.0 W 
2,635 

UT4 Yepun 8.2m 

NACO/H 

0 .2/0.2 

G. Hau 

San Pedro 

de Atacama 

22:57:12.3 S 

68:10:47.6 W 
2,397 

Caisey 0.5m f/8 

CCD/V 

3/4.58 

A. Maury 



Caisey 0.5m f/6.8 

CCD/B 

4/4.905 

L. Nagy 



CAO 4 0.4m 

CCD/R 

4/6.35 

J.F. Soulier 


— 

ASH2 0.4m 

STL11000 
technical problem 

N. Morales 



OPSPA 5 0.3m 

CCD/clear 

5/11.1 

J. Fabrega Polleri 

Huancayo 

12:02:32.2 S 

75:19:14.7 W 
3,344 

Meade 8 0.2 m 

CCD/clear 

10.24/10.24 

Negative chord 

E. Meza 


Notes. 1 PROMPT: Panchromatic Robotic Optical Monitoring and Polarimetry Telescopes. 

. 3 B&C: Boiler & Chivens. 4 CAO: Campo Catino Observatory. 
5 OPSPA: Observatorio Panameno en San Pedro de Atacama. 
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18 July 2012 



Time 


Fig. 3.— The blue curves are a simultaneous fit to the 18 July 2012 light-curves, using the best atmospheric 
model described in Table [4] and Fig. 13 The number at the lower right of each panel is the value of 
(Eq. A8), i.e. the x 2 per degree of freedom for each corresponding fits. Each panel spans 3 minutes of data, 
with the vertical tick marks located at 04:13 UT. All the light curves show the total flux (star+Pluto+Charon) 
plotted at the same vertical scale. The horizontal bars on the Cerro Burek, San Pedro de Atacama and 
Huancayo data points represent the respective integration times. The zero flux is indicated by the solid 
horizontal line at the bottom of each panel, together with the residuals (data minus model). The dotted 
horizontal lines mark the fitted zero stellar fluxes (or equivalently, the Pluto+Charon contribution to the 
total flux), obtained using our best Pluto atmospheric model. The blue horizontal line in the Paranal panel 
marks the measured zero stellar flux at that station, the only one at which a photometric calibration was 
possible (see text and Figs. In the Paranal panel, we have also added the residuals (labelled “iso.”) 

obtained by forcing an isothermal mesosphere at T- lso = 95.5 K. The residuals have been averaged over 5-s 
time intervals and shifted vertically by -0.12 for better showing the clear discrepancy between the isothermal 
mesospheric model and the data. Other values chosen for Tj so would result in the same qualitative behavior. 
In essence, isothermal mesospheres do not provide satisfactory fits to the NACO light-curve. 
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04 May 2013 



Time 


Fig. 4.— The same as Fig. [3] for the 04 May 2013 event. Each panel now represents 6 minutes of data, 
with the vertical tick mark located at 08:22 UT. Note that the two light-curves from San Pedro ( “SP”) de 
Atacama have been displaced vertically by ±0.1 for better viewing. 
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04 May 2013 



Fig. 5.-- The same as Fig. [4j but for the light-curves that were not included in the fit, either due to lower 
SNR, or interruptions during the acquisition. See Table [5] for intrumental details (“SP” refers to San Pedro 
de Atacama and acronyms refer to telescope used in that station.). Note that the Leoncito, Danish and 
SP light-curves duplicate the observations of the Cerro Burek, La Silla TRAPPIST and Caisey telescopes, 
respectively. 
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Fig. 6.— Photometric calibration of the July 18, 2012 event (Paranal/VLT, NACO H-band). Left - Image 
taken some 20 minutes before the event, showing the small separation between Pluto, Charon and the star 
(~ 1”). Right - The same image after a digital coronagraphy treatment that removed the stellar image. See 
text for details. 
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18 July 2012 NACO 



d ioo 150 200 250 
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Fig. 7. — Left - details of the fit to the NACO 18 July 2012 light-curve (see also the middle-left panel of 
Fig# The horizontal blue line in the gray shaded area indicates the Pluto + Charon contribution to the 
total observed flux and its 1 -a error bar, 0.1184 ±0.007. The residuals curve at the bottom clearly shows the 
spike activity at ingress and egress. Right - Expanded view of the left panel. The data have been binned over 
1 second-time intervals to better show the flux decrease during the central phase of the occultation. The flux 
of the primary stellar image is plotted in green, while the blue curve is the sum of the primary and secondary 
images, according to the model (see Fig. [2] and Appendix). Thus, the contribution of the secondary image 
is the difference between the blue and green curves. Note the interruption of data acquisition (about 3 s) 
at mid-occultation, necessary to the writing of the data cube before the start of the next data cube. The 
inclined gray line is a linear fit to the central part of the light-curve, which illustrates the ingress/egress 
asymmetry of the residual stellar flux. The vertical axis inside the box at left indicates the value of the 
residual, normalized stellar flux. It shows that the stellar flux decreased from about 2.3% to about 1.8% of 
its full unocculted value during that interval. The systematic error on those values is ±0.8% (corresponding 
to the shaded area). 
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Fig. 8.— In all the panels, the black solid line is the model that best fits all the 18 July 2012 NACO 
light-curves, see Figs. [3] and Table |4j The red (resp. blue) lines are the particular profiles obtained from the 
inversion of the NACO 18 July 2012 light-curve at ingress/summer (resp. egress/winter). The shaded areas 
at the top of the profiles indicate the expected ±ler fluctuations caused by the photometric noise, see text. 
The shaded areas at the bottom of the profiles are the ±1 <t uncertainty domain caused by the uncertainty 
on the Pluto + Charon contribution to the 18 July 2012 NACO light-curve, see Fig. [7] (a) Molecular density 
vs. radius (assuming a pure N 2 atmosphere); (b) temperature vs. pressure; (c) temperature vs. radius; (d) 
temperature gradient vs. radius. The two gray temperature profiles in panels (b), (c) and (d) show the effect 
of different temperature boundary conditions for the egress NACO profile. More precisely, those profiles 
differ from the nominal one (blue lines) by ±5 K at 1,390 km. The oblique solid line at the left of panel (b) 
is the vapor pressure equilibrium limit for N 2 (Fray and Schmi tt|20 09). Nitrogen should condense at the left 
of this line, so that its intersection with the temperature profile may define Pluto’s surface in the absence of 
troposphere, see Fig. [9] and text for details. The dash-dotted line in panel (d) is the dry adiabat for a pure 
N 2 atmosphere. 
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T (K) 


dT/dr (K km 1 ) 


Fig. 9.— Left - Expanded view of the bottom of the temperature profiles shown in Fig. [8] The bullet is 
the intersection with the nitrogen condensation line. The error bar attached to its positions is defined by 
the radial extension of the shaded uncertainty domain. Right - The corresponding expanded view for the 
temperature gradient. 
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Fig. 10. The temperature profiles derived from the inversion of our best occultation light-curves obtained 
on 18 July 2012 and 04 May 2013. The dotted line is our global, best-fitting temperature profile (also shown 
in Figs. [S] and 


131 . 
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Fig. 11. - Calculation of cooling rates by CO and HCN. The y-axis is the distance from Pluto’s center, with 
the surface position assumed here at 1,184 km. Left panel - Cooling rates assuming the thermal profile from 
this work. Red and green curves: CO cooling rates for qco = 5 x 10 -4 and 200 x 10 -4 , respectively. The 
other three colored curves show the HCN cooling rate for the corresponding HCN profiles. Right panel - 
HCN mixing ratios profiles. The black and purple curves make use of the thermal profile from this work. 
Due to the significantly cold temperatures above ^1300 km, HCN is limited by saturation throughout the 
atmosphere, except in a limited region at 1210-1270 km for an assumed ^hcn = 10 -7 . The blue curve shows 
the hypothetical case of a uniform (i.e. non-limited by saturation) 5 x 10 -5 uniform HCN mixing ratio. 
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t (k) M'g'/Mg 


Fig. 12.— Left panel - Solid line: our best temperature profile (see Fig. 13). Gray line: an example of a 
prescribed profile with milder mesospheric thermal gradient, here a 10 K drop between the stratopause and 
the mesopause. Right panel - The ratio g! g' /ng, as defined by Eq. [2j corresponding to the gray, prescribed 
profile of the left panel. The points numbered 2, 3 and 4 correspond respectively to the stratopause, the 
inflexion point and the mesopause (see also Fig. 13 and text for details). 
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A. Synthetic light curves 
A.l. Parametrized temperature profile 


We define a parametric model for Pluto’s atmosphere temperature profile, T(r), where r is the radius, 
i.e. the distance to Pluto’s center. The model must be detailed enough to capture the main features revealed 


by the inversions (Figs. 10 and 13), but still simple enough to allow an easy and meaningful control of T(r) 


and an assessment of the error bars associated with each parameter. The features we want to describe are: 
(?) a thin stratosphere just above the surface, with a strong increase of temperature with altitude, ( ii ) an 
“elbow” where the temperature reaches its maximum, marking the stratopause, (in) an intermediate region 
with a mild negative gradient, and finally ( iv ) an isothermal upper branch. 

These features define three regions, bounded by four points 1,...4 at prescribed radii r\, ...n, see Fig.[l3| 
More precisely, the profile T(r) is generated as follows: 


cl ■ r + c2 ■ T(r) + c3 • rT(r ) + c4 • r 2 + c5 • T 2 (r) = 1 
< T(r) = c6 + c7 ■ r + c8 ■ r 2 + c9 ■ r 3 
. T(r) = T iso 


for rq < r < r^ (hyperbolic branch) 
for C 3 < r < T 4 (polynomial branch) 
for C 4 < r (straight line) 


(Al) 


Note that r 2 does not appear in the equations above, and is defined as the radius where the temperature 
reaches its maximum (Fig. 13). The functional forms chosen here (hyperbolic, polynomial and straight lines) 
are not based on physical grounds, but rather, are empirical and simple formulae that satisfactorily fits the 
observed profiles (Fig. [ 8 ]). 

The parameters cl, ...c9 are determined to ensure that T(r ) is continuous both in temperature and its 
derivative, dT/dr , at points 1, 3 and 4. Those conditions provide algebraic systems that are solved by a 
classical Gauss-Jordan method (Press et al. 19921. 


In practice: (1) we fix the temperature Ti = T(ri) at the bottom of the profile, together with its 
gradient ( dT/dr) i . (2) We fix the value of the maximum of temperature T 2 = T(r 2 ) at r 2 and the temperature 
T 3 = T(rs) at the inflexion point 3. We thus have three boundary conditions for T: Ti, T 2 , T 3 at ri, r 2 and r 3 , 
respectively, and two boundary conditions for dT/dr: (dT/dr)(n) = (dT/dr) 1 and (dT/dr) (r 2 ) = 0, which 
fixes the five coefficients cl,...c5. Note in passing that the values of cl, ...c5 then impose the temperature 
gradient (dT/dr )3 at r^,\ (3) We fix the temperature T- lso at r±, the point where the isothermal branch 
is reached. This provides two boundary conditions in T: T 3 and T lso at r$ and r±, respectively, plus 
two boundary conditions for dT/dr : (dT/dr )3 at r% and (dT/dr)(r^) = 0, thus fixing the remaining four 
coefficients c6...c9. 


The locations of points 1, 2, 3 and 4 in the space (T,r) are chosen to best fit the observed profiles, see 
the main text for details. Once T(r) is defined, the gas number density profile n(r ) is obtained by integrating 
the first order differential equation: 

1 dn 
n dr 


d9(r) 

kT 


1 

T 


dT 

dr 


(A2) 


derived from the equation of state for an ideal gas, and the hydrostatic equation. Here, 


g(r) = 


GM 

r 2 


(A3) 
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Table 3: Parameters of the two occultations 



18 July 2012 

04 May 2013 

Star coordinates 1 

a= 18h 32m 14.6720s 

5= -19d 24’ 19.295” 

a= 18h 47m 52.5322s 

-19d 41’ 24.3738” 

Ephemeris 

DE413/PLU022 

DE413/PLU031 

Pluto geocentric distance 

4.68244 x 10 9 km (at 04:13 UT) 

4.76882 x 10 9 km (at 08:23 UT) 

Sub-observer and sub-solar latitudes 2 

B= +47.10d, B’= +47.54d 

B= +49.95d, B’= +48.64d 

Pluto’s north pole position angle 2 

P= -56.88d 

P= -52.91d 

Shadow velocity 

ss 23.0 km s -1 

ss 10.6 km s -1 

Magnitudes 3 

V= 14.7, R=13.7, K= 10.9 

V= 14.1, R=14.0, K= 12.4 


1 J2000, UCAC2 system. 

2 Assuming the Pluto’s north pole position (J2000) of Tholen et al. (2008): a p — 08h 52m 12.94s, <5 P = -06d 

10’ 04.8” _ 

3 From NOMAD catalog (Zachar i'as et al.| |2004) 


Table 4: Pluto atmospheric model 


Pluto’s mass 1 
Nitrogen molecular mass 2 
Nitrogen molecular refractivity 3 
Boltzmann constant 


Physical parameters 
GM = 8.703 x 10 11 m 3 s -2 
/i = 4.652 x 10 -26 kg 

K = 1.091 x 10 -23 + (6.282 x 10 -26 /A 2 m ) cm 3 molecule -1 
k = 1.380626 x 10 -23 J K -1 


r l, Ti, dT/dr(ri) 
r2, T 2 
T3, T 3 

Tj, T 4 _ 

cl, c2 
c3, c4 
c5, c6 
c7, c8 
c9 


The nine free parameters of the best temperature profile 4 
1,190.4 ± 1 km, 36 K, 16.9 K km -1 

1.217.3 km, 109.7 K 

1.302.4 km, 95.5 K (implying dT/dr(r 3 ) = —0.206 K km -1 ) 
1,392.0 km, 80.6 K 

1.41397736 x 10 -3 , 2.59861886 x 10 -3 
-2.19756021 x 10 -6 , -4.81764971 x 10 -7 
8.66619700 x 10 -8 , -3.6213609 x 10 4 
8.2775269 x 10 1 , -6.27372563 x 10 -2 
1.58068760 x 10 -5 


The three free parameters particular to each event 5 


18 July 2012 

04 May 2013 

Pressure at r = 1,275 km, pi ,275 

2.16 ± 0.02 /ibar 

2.30 ± 0.01 /xbar 

Time of closest geocentric approach 

04:13:37.24±0.07 UT 

08:22:27.11±0.09 UT 

Distance of closest geocentric approach 6 

-404.6 ± 2.7 km 

-723.5 ± 2.7 km 


^Tho len et al. (2008). 2 Assumed to be the only constituent in the ray tracing code, see text. "Washburn 
(1930). For both NACO observations of 2012 and 2103, H band (A = 1.6 /j,m) was used. 4 Or equivalently, 
the nine coefficients cl, ...c9, see text and Fig. 13 for the definition of the various quantities given here. 5 So 
that there are twelve free parameters for each date. 6 Distance of Pluto’s center to the star at closest 
approach, projected in the sky plane, as seen from the geocenter. A negative value means that Pluto’s 
center went south of the star in the sky plane. 
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is the acceleration of gravity, assuming a spherical, homogeneous planet. The values of /i (Table [4]) corre¬ 
sponds to molecular nitrogen, assumed to be the unique gas present in the atmosphere. Also listed in Table|4] 
are the Boltzmann constant k and Pluto’s mass parameter GM. 


A boundary condition is required to integrate Eq. X2 e.g. the pressure £> 1,275 at r = 1,275 km, which 


fixes the needed boundary condition 711,275 = Pi, 275 /^^ 1 , 275 . Finally, the refractivity v{r) of the gas (index 
of refraction minus unity) is given by 

v{r) = K ■ n(r ), (A4) 


where the molecular refracticity is given in Table [3j assuming again pure molecular nitrogen. Once v{r) 
is obtained, we can derive the vertical refractivity gradient dv/dr that is used in the ray tracing code, see 
below. 


The inversions proceed the other way around: the light-curves provide dv/dr{r ) through an abelian 
integral (Vapillon et al. 1973), then v{r), from which n(r) is derived (Eq. A4), followed by the temperature 
profile, once a boundary condition is given for T (Eq. A2) 


A.2. Ray tracing 


For small values of v (as it is the case here) and under spherical symmetry assumption, a stellar ray is 
deviated by du> = {dv/dr) ■ ds (Snell-Descartes law) as it moves along an elementay path ds. In principle, a 
ray tracing code should account for the curvature of the stellar ray as it is refracted in the atmosphere. In 
practice, however, it is enough to assume that the ray has a rectilinear trajectory in the entire atmosphere. 
In fact, the maximum total deviation ui suffered by the ray is very small for ground-based occultations, more 
precisely of the order of Pluto’s apparent angular radius, ~0.05 arcsec, so that w <~ 3 x 1CT 7 rad. Most of 
that deviation o ccurs in the deepest sca le height H traversed at radius r, which represents a traveled length 
of l ~ -\J2tt rH (Baum and Code 1953). Taking typical values of r ~ 1,200 km and H <~ 50 km, we get 
l <~ 600 km, i.e. a deviation inside the atmosphere of ~ oj ■ l < 0.2 meters, which is negligible compared to 
the scales probed by ground-based stellar occultations. 


The numerical integration of Eq. A2 using a second order scheme, provides ro(r;) at discreet layers of 
radii ?y, from which the refractivity Vi and its gradient {dv / dr){r/) are calculated. The total deviation along 
the straight line s is then: 

to = A u>i = ^^{dv/dr){ri) ■ As,, (A5) 

i i 

where As^ is the path along s traveled inside the layer i. Then, for a closest approach r of a ray to Pluto’s 
center, the corresponding distance 2 to the shadow center upon arrival on Earth is 


z = r + oj{r) ■ D , 

where D is Pluto’s geocentric distance. The observed stellar flux is then 

= 4 = 1 1 


+ Ddoj/dr ’ 


(A 6 ) 


(A7) 


where f = r/z is the focusing factor due to the (assumed circular) limb curvature, see Sicardy et al. (1999). 


The thickness A?y of the individual refracting layers has been adjusted to 30 meters to minimize numer¬ 
ical noise, while keeping computing times reasonably low. Similarly, the sampling for r (the closest distance 

















Radius (km) 
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of the rays to Pluto’s center) has been adjusted so that adjacent rays arrive at separation A z « 1 km in the 
shadow plane. 


Once the table (r,z,$(z)) has been completed, the synthetic flux at a given site and given moment 
(corresponding to a distance z 0 bs °f the observer to the shadow center) is calculated by interpolation. If 
several stellar images are present, all the fluxes are summed. In the particular case of a spherically symmetric 
atmosphere, and for a given distance z 0 bs, there is a primary image corresponding to z = z 0 bs> and a secondary 
image corresponding to z = —z 0 b s . 


The lowest radius n considered in the model (1,190.4 km, see Tableland Fig. 13) is adjusted so that the 
corresponding flux received in the shadow plane is ~ 10 -3 of the unocculted stellar flux, negligible compared 
to the noise level of the best light-curves. The upper limit for the atmosphere has been fixed to a radius of 
about 2,300 km. This corresponds to a pressure level of about 0.05 nbar, at which point the stellar drop is 
several orders of magnitudes less than the noise in our best light-curves. 




Temperature (K) z (km) 


Fig. 13.— Left - The temperature profile T{r) that best fits our 2012 and 2013 light-curves, see Figs. [3] and 
|4j The parameters used to generate this profile are given in Table [4] Total thickness of the atmosphere: 
1100 km, vertical sampling of the model: 0.03 km. Right - The corresponding synthetic flux in the shadow 
plane for the 18 July 2012 occultation. Here, z is the distance to the shadow center, with the four points 
corresponding to those of the left panel. 


The profile that best fits our light-curves is shown in Fig. |13| The trajectories of the primary and 
secondary stellar images as seen from VLT on 18 July 2012 are displayed in Fig. [2] 

The best fits to the observed light-curves are shown in Figs. [3] and [4j Their quality is assessed through 
the so-called y 2 per degree of freedom: 


,2 _ X _ y 

Xd°f — N _ M - N _ M 2-^ 


1 


N 


N-M 


i=1 


Tub-.: T- 


(A8) 


where 4> obSj i (resp. <I> syn ,i) is the observed (resp. synthetic) stellar flux of the i th data point, cq is its associated 
standard deviation, where the summation is extended to the N data points from all the light-curves used in 
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the fit, and M is the number of free parameters of the model. As we have nine coefficients Ci, ...eg to define 
T(r) (Eg. |A1[), a boundary condition pi ,275 and two quantities to define Pluto’s center, M = 12. 


B. Noise propagation 
B.l. Photometric noise 


Here we estimate the effect of photometric noise in an occultation light-curve on the retrieved density, 
temperature and temperature gradient profiles. We denote <5 the fluctuation of a given quantity, and a = Vs 2 
its standard deviation, where the bar denotes average values. For estimation purposes, it is enough to assume 
here (but not in the ray-tracing or inversion procedures) that the atmosphere has locally a constant density 
scale height H that is small compared to the planet radius. The stellar flux is then given by the Baum and 
Code (BC) equation (Baum and Code||1953): 




z 1/2 


H 


(Bl) 


where Z\/ 2 is the position in the shadow plane at which $ = 1/2 (half-light level). 

We focus on the top of the profiles, corresponding to $ ~ 1, so that Eq. Bl becomes $( 0 ) ~ 1 — 
exp[— (z — Z 1 / 2 )/ H], Morever, for <f> ~ 1, the stellar ray deviation u> is small, and we can equate r and z, 
see Eqs. 


A6 


and 


A7 


where / ~ 1. In the BC approximation, we have u> ~ —v^2nr/H, where v is the 
refractivity at r. As the atmosphere density profile is basically exponential, duj/dr ~ —c o/H = u^/2nr/H 3 , 


so that Eq. |A7| can be used to estimate the expected refractivity corresponding to a stellar flux d>, namely 
v ~ (1 - <S>)^H 3 /2irrD 2 

We denote v Q and r 0 the refractivity and corresponding radius where the stellar drop is equal to the 
standard deviation of the flux, < 7 $, i.e. 


z'o ~ 0$ 


H 3 

2iTrD 2 


(B2) 


Thus, ro is the radius where the stellar drop starts to be barely significant, given the photometric 
noise. At the upper part of the profiles, we have H ~ 60 km. The 18 July 2012 NACO lightcurve has 
a photometric standard deviation of er$ = 0.011. Using the value of D given in Table [3j we obtain v 0 ~ 
1.3x 10 -11 . Assuming a pure N 2 atmosphere, we obtained the corresponding molecular density no = ~ 

6 x 10 12 cm -3 , which is reached at radius 


ro ~ 1, 565 km. 


For <!> ~ 1 (and / ~ 1), and using the results above, Eq. A7 provides $ ~ 1 —Ddu/dr = l+y / 27r rD 2 /H(dv/dr). 


For a noise-free light-curve, we expect $ = 1 — cr$ exp[—(r — rf)/H], In reality, d* is affected by fluctuations 
so that the retrieved refractivity gradient is in fact: dv/dr = (y 0 /H) ■ [—exp[—(r — r 0 )/H\ + 8&/(j\. 
Consequenly, the standard deviation associated with each point of the ( du/dr)(r ) profile (and restricting 
ourselves to the top of the profile) is: 


@dvI dr 


vo 

H 


The profile ( du/dr)(r ) is the primary result derived from the light-curve, and from which all the 
other profiles are deduced. Once dv/dr is known, we have to estimate v{r) = V\ + J^(dv/dr)dr, where 
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(t’i, = v{r\)) is a boundary condition. The integration is performed numerically by taking v(r) = 

v(r\) + ’Y/i—\{du/dr){ri) ' Ar, where Ar is the spatial sampling of the data (i.e. A r = the star velocity 
perpendicular to the limb multiplied by the exposure time). Thus, r* = r\ + (i — l)Ar and N = \r — ri|/Ar. 
Adding the variances associated with individual (di//dr)(7y)’s, we obtain: 



where the second approximation stems from the fact that rq is chosen close to Tq and that we are considering 
here the few top scale heights of the profiles, so that |r — ri\/H ~ 0(1). Note that a v = 0 for r = rq. This 
is because (rq , iq) is an arbitrary boundary condition, and such, has no associated error bars. 

From n = v/K, we obtain the standard deviation associated with the density gradient and the density 
itself: (jdn/dr = cr dv/dr/K and cr n = a v /K. Moreover, from Eq. A2 and assuming an isothermal upper 
atmosphere, we obtain S^T/dr = ~{T/n)8d n /dr: so that (JdT/dr = (' T/n)ad n /dr■ Finally, the temperature 
profile is obtained from the numerical integration of T(r) = T 2 + J^(dT/dr)dr, where (= T(r 2 )) is 
an arbitrary boundary condition. Using the same line of reasoning as for n(r), we obtain rr t by adding 
the variances &d T / dr of all the points i = 1...N involved in the integration, where now N = \r — 7'2|/Ar. 
Combining the results above, we obtain the following standard deviations for n(r), T(r) and ( dT/dr)(r ): 



(r2-ro)/H 


a dT/dr 


H C 


( r-r 0 )/H 


(B3) 


Fig-! shows the ±lcr envelopes at the upper parts of the various profiles. We take here r 2 = 1, 390 km, 
the radius at which we fix a prescribed temperature T 2 ~ 81 K. Note again that or = 0 at r = V 2 , as (r 2 , T 2 ) 
is an arbitrary boundary condition. Finally, the envelopes ±lcr are plotted only down to the half-light level 
(r ~ 1, 290 km), as the estimations made here apply only for the upper part of the light-curve. In any case, 
below that level, the uncertainties in the profiles are dominated by the uncertainty on the background Pluto 
+ Charon contribution, see below. 


B.2. Effect of the Pluto and Charon flux contributions 


The stellar flux reaches its minimum value in the shadow at typically z m - m ~ (^i/ 2 )/2, i.e. half-way 
between the half-light level and the shadow center, where the central flash occurs (Fig. [Tdj). At the minimum, 
we have from Eq. B1 $ m i n ~ H/(zi /2 — ~ m in) ~ 2H/z\/2- Eq. A2 then provides 


H = 


dn/dr 


T 

(iJ-g/k) + ( dT/dr) 


z l/2 


d> 

^ n 


(B4) 


At the bottom of the temperature profile (stratosphere), i-ig/k and dT/dr are of same order of magnitude. 
Consequently, increasing the value of the Pluto + Charon contribution to the light-curve decreases the value 
of >!> (Fig. [7]), thus increasing the retrieved gradient dT/dr. This is illustrated in Fig. [9j 


B.3. Effect of initial conditions 


Once the density profile n{r) is derived from the inversion, Eq. A2 yields the temperature profile T(r), 
provided a boundary condition X), = T(rf,) is fixed at an arbitrary level 7v Let us consider two possible 
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solutions T(r) and (T+AT)(r) that differ by A T(rb) at 7 - 5 , then expanding Eq. A2 to first order in (AT/T)(r), 
we obtain: 


d 

(AT\ 

~ 1 1 

( at\ 

dr 1 

\ T )' 

H ' 

l T ) 


where we have approximated H ~ kT/gg. Thus, as r increases, the relative difference A T/T diverge 
exponentially as: 


( r ) „ (** 


V T 


\ T 


in) ■ e 


(' r-r h )/H 


(B5) 


This exponential divergence should not been confused with the one that is provided by Eq. B3 for <jt- 
The latter tends to zero as the noise tends to zero, while the former is inherent to the nature of Eq. |A2[ 




